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Abstract 

We  consider  parameter  estimation  in  linear  models 
when  some  of  the  parameters  are  known  to  be  integers. 
Such  problems  arise,  for  example,  in  positioning  using 
phase  measurements  in  the  global  positioning  system 
(GPS.)  Given  a  linear  model,  we  address  two  problems: 

1.  The  problem  of  estimating  the  parameters. 

2.  The  problem  of  verifying  the  parameter  esti¬ 
mates. 

Under  Gaussian  measurement  noise: 

•  Maximum  likelihood  estimates  of  the  parameters 
are  given  by  solving  an  integer  least-squares  prob¬ 
lem.  Theoretically,  this  problem  is  very  difficult 
to  solve  {^f'P~haxd.) 

•  Verifying  the  parameter  estimates  (computing 
the  probability  of  correct  integer  parameter  es¬ 
timation)  is  related  to  computing  the  integral  of 
a  Gaussian  PDF  over  the  Voronoi  cell  of  a  lattice. 
This  problem  is  also  very  difficult  computation¬ 
ally. 

However,  by  using  a  polynomial- time  algorithm  due  to 
Lenstra,  Lenstra,  and  Lovasz  (LLL  algorithm): 

•  The  integer  least-squares  problem  associated 
with  estimating  the  parameters  can  be  solved  ef¬ 
ficiently  in  practice. 

•  Sharp  upper  and  lower  bounds  can  be  found  on 
the  probability  of  correct  integer  parameter  esti¬ 
mation. 

We  conclude  the  paper  with  simulation  results  that  are 
based  on  a  GPS  setup. 

Key  words:  Integer  parameter  estimation,  Linear 
model.  Integer  least- squares,  GPS. 
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1  Problem  Statement 

1.1  Setup 

Throughout  this  paper,  we  assume  that  the  observation 
y  e  R”  is  related  to  the  unknown  vectors  x  {real) 
and  z  {integer)  through 

y  =  Ax  Bz  -\-v,  (1) 

where  A  €  (full  column  rank)  and  B  G  R”^^^ 

are  known  matrices.  The  measurement  noise  i;  G  R'^  is 
assumed  to  be  Gaussian  with  zero-mean  and  covari¬ 
ance  matrix  identity,  i.e.,  Ar(0,/nxn)-  This  can  be 
assumed  without  loss  of  generality,  as  we  can  always 
rescale  equation  (1)  by  the  square  root  of  the  covari¬ 
ance  matrix  of  v.  Given  y,  A,  and  R,  our  goal  is  to  find 
and  verify  estimates  of  the  unknown  parameters  x  and 
z.  Some  previous  results  are  given  in  [6]  and  [7]  and 
references  therein. 

1.2  Estimation  problem 

We  consider  maximum  likelihood  (ML)  estimates  x^i. 
and  Zml  for  x  and  z  respectively  that  maximize  the 
probability  of  observing  y,  z.e., 

{xml,Zml)=  argmax  {pY\x,z{y\x,z)] .  (2) 
(x,z)  eRPxZ? 

1.3  Verification  problem 

Since  z  is  an  integer-valued  vector,  we  have  the  chance 
of  estimating  (or  detecting)  it  exactly.  As  a  result, 
for  the  verification  step,  we  compute  the  probability  of 
estimating  z  correctly,  i.e., 

Pe  =  Prob(zML  =  ^),  (3) 

Clearly,  the  larger  this  value  the  higher  the  reliabil¬ 
ity  on  the  estimates  Xml  and  Zml-  Another  reliability 
measure  can  be 

Prob(||x -XmlII  <  e),  (4) 

where  ||  •  ||  denotes  the  2-norm.  This  gives  us  the  prob¬ 
ability  of  the  real  parameter  estimate  lying  in  a  ball 


of  radius  e  of  its  actual  value.  A  combined  reliability 
measure  is 

Prob(||a:  -  rcMLlI  <  ^  &  2;ml  =  2)  = 

Prob(j|a:  -  Xmi.\\  <  e\  z^l  =  z) »  Prob(zML  =  z). 

(5) 

In  practical  cases,  e  is  small  only  if  Zml  =  -2:,  and  there¬ 
fore,  (4)  and  (5)  are  almost  equal. 


2  Problem  formulation 

Since  v  is  Gaussian  we  get 


Also,  (8)  can  be  written  in  terms  of  G  and  y  in  the 
following  equivalent  form 

Zuv  =  argmin  ||^  -  Gzf.  (12) 

zeZ^ 

The  set  {Gz  |  2:  €  Z*}  is  a  lattice  in  R*.  Equation  (12) 
states  that  z^l  is  found  by  computing  the  closest  lattice 
point  to  the  vector  y.  In  addition,  according  to  (11), 
y  is  off  from  Gz  by  u.  Therefore,  as  long  cis  u  is  small 
enough  such  that  y  remains  closer  to  Gz  than  anyother 
point  in  the  lattice,  the  estimate  of  2:  is  correct  (i.e., 
Zml  =  z.)  This  is  equivalent  to  remaining  inside 

the  Voronoi  cell  of  the  lattice  point  Gz'^,  Thus 


(xml,  ^ml)  =  argmin  Hy  -Ax-  Bzf, 
(x.z)  €RPxZ« 


Using  completion  of  squares  arguments  it  can  be  shown 
that 


(a:ML,2ML)=  argmin  |(a:  -  ^  (i  -  Xj^,)  + 

(x,^)  e  BPxZ" 

(z  -  z)^T“^(z  -  z)  +  y^Ay}  , 

(6) 

where 


T 

’|2 

z 


{A'^Ar\ 

[I  -  AEA^)  B)  \ 
EA^iy  -  Bz), 

TB'^  {I -AEA^)y, 


Pc  =  Prob(M  €  Vb)  where  u  ^  N{Q,Igxq),  (13) 

and  Vo  is  the  Voronoi  cell  of  the  origin.  According 
to  (10),  X|j  =  x  +  u)  where  w  ~  iV’(0,H),  and  if  z  =  Zml 
we  have 

Prob(||x  -  XmlII  <  e  1 2  =  Zml)  =  Prob(||E^/^wll  <  e), 

(14) 

which  is  equal  to  the  probability  of  w  falling  inside  the 
ellipsoid  <  e. 

Formulas  (8),  (9),  (12),  (13),  and  (14)  give  the  esti¬ 
mates  and  their  measures  of  reliability. 


3  Verification  Problem 


and  A  >  0  is  a  constant  matrix.  Clearly,  from  (6) 

2ml  =  argmin  {{z  -  z)'^i:~^iz  -  z)}  ,  (8) 

zeZ^ 

^ML  —  ^|2mL* 

As  a  result,  2:ml  (and  therefore  Xml)  is  given  by  the 
solution  to  a  least-squares  problem  involving  integer 
variables.  Because  of  the  integer  nature  of  z,  comput¬ 
ing  Zmi  is  very  difficult  and  is  known  to  be  VP-hard 
(cf.  [IJO 

In  order  to  calculate  the  measures  of  reliability  Pc  and 
Prob(l|x  -  x||  <  e|2;  =  Zml)  we  need  to  know  the 
probability  distributions  of  our  estimates.  Note  that 
we  can  easily  read  the  covariance  matrices  of  z  and  x\z 
from  (6)  as 

cov(i)  —  T  and  cov  (x|^)  =  5. 

Therefore 

i~A'(z,T)  and  x\z^N{x,E)j  (10) 

so  for  example  z  =  z-\ru  where  u  ~  ^(0,  T).  Multiply¬ 
ing  both  sides  by  G  :=  (the  whitening  matrix  of 

u)  and  by  defining  y  :=  Gz  we  get 

y=^Gz-\-u^  where  u  ^  N{Q,Iqxq)’  (H) 


As  noted  in  §1,  the  probability  of  detecting  the  integer 
parameter  z  correctly.  Pc  in  (13),  can  be  used  to  verify 
the  ML  estimates.  The  reliability  measure  given  in  (5) 
also  requires  the  calculation  of  Pc.  Therefore,  an  es¬ 
sential  step  in  verifying  Xml  and  Zml  is  the  calculation 
of  Pc. 

It  turns  out  that  computing  Pc  is  very  difficult  compu¬ 
tationally  and  therefore  we  are  motivated  to  find  fast 
ways  to  approximate  Pc  that  enable  real-time  reliability 
tests  instead.  In  fact,  as  we  will  see,  easily  computable 
upper  and  lower  bounds  on  Pc  exist  which  become  tight 
as  Pc  approaches  unity. 

In  §3.1  we  see  that  the  choice  of  G  for  a  given  lattice 
is  highly  nonunique.  Although  Pc  only  depends  on  the 
lattice  and  not  on  the  specific  choice  of  G,  this  is  not 
true  for  the  upper  and  lower  bounds  on  Pc.  There¬ 
fore,  we  can  sharpen  these  bounds  by  optimizing  over 
the  family  of  admissible  matrix  representations  for  the 
lattice.  We  introduce  the  LLL  algorithm  that  is  an  ex¬ 
tremely  useful  tool  for  sharpening  bounds  on  Pc,  and  as 
will  become  clear  later,  it  is  also  useful  for  improving 
the  efficiency  of  the  algorithm  for  solving  the  integer 
least-squares  problem  for  computing  Zml- 

^The  Voronoi  cell  of  a  point  Gz  in  a  lattice,  is  the  set  of  all 
points  in  space  that  are  closer  to  Gz  than  anyother  point  in  the 
lattice. 


3.1  Lattice  Bases 

Let  us  denote  a  lattice  L  G  generated  by  the  matrix 

G  =  [fli52---5,]eR«^’by 

L  =  L(G)  =  {Gz\ze  Z«}  . 

The  set  of  vectors  {pi,  <9^2?  •  •  •  >  Pg}  is  called  ‘a’  basis  for  L 
since  all  lattice  points  are  linear  combinations  of  integer 
multiples  of  these  vectors.  We  say  ‘a’  and  not  ‘the’  be¬ 
cause  a  lattice  L  has  infinitely  many  bases.  The  lattice 
generated  by  the  matrix  GF,  where  F  is  any  unimodu- 
lar  matrix,  and  that  of  G  are  identical.  Such  matrices 
F  have  the  property  that  the  elements  of  F  and  F~^ 
are  integers  (or  equivalently,  the  elements  of  F  are  in¬ 
tegers  and  detf  =  ±1.) 

Of  the  many  possible  bases  of  a  lattice,  one  would  like 
to  select  one  that  is  in  some  sense  nice  or  simple,  or 
reduced.  There  are  many  different  notions  of  reduced 
bases  in  the  literature;  one  of  them  consists  of  requir¬ 
ing  that  all  its  vectors  be  short  in  the  sense  that  the 
product  IIpiII  ||p2||  *  •  •  ||p^||  is  minimum.  This  problem 
is  known  as  the  minimum  basis  problem  and  has  been 
proved  to  be  J\fV‘haid  (cf.  [1].) 

Although  the  minimum  basis  problem  is  jV^'P-hard  there 
exists  a  polynomial-time  algorithm  for  computing  a  ba¬ 
sis  for  a  given  lattice  that  is  in  some  sense  reduced. 
This  algorithm  is  due  to  A.  K.  Lenstra,  H.  W.  Lenstra, 
Jr.,  and  L.  Lovasz  (LLL)  and  is  very  efficient  in  prac¬ 
tice.  Given  a  lattice  generator  matrix  G,  the  LLL  al¬ 
gorithm  gives  a  new  lattice  generator  matrix  GF  with 
F  unimodular  such  that  GF  has  two  nice  properties. 
Roughly  speaking  these  are: 

•  The  columns  of  GF  (the  new  basis  vectors)  are 
almost  orthogonal. 

•  The  2-norm  of  the  columns  of  GF  (the  length  of 
the  new  basis  vectors)  are  not  arbitrarily  large. 

For  an  exact  discussion  of  these  properties  refer  to 
Chapter  5  of  [1].  The  LLL  algorithm  is  described  in 
the  appendix. 

3.2  Calculating  Pc 

Computing  Pc  as  given  in  (13)  requires  the  calculation 
of  a  suitable  description  of  Vo  the  Voronoi  cell  of  the 
origin  in  the  lattice  generated  by  G.  This  is  possible  of 
course,  but  it  is  very  difficult  computationally  (cf.  [8].) 
Adding  to  this  the  computational  cost  of  the  numeri¬ 
cal  algorithm  to  perform  the  integral  of  the  Gaussian 
probability  density  function  over  Vo,  we  conclude  that, 
although  the  exact  calculation  of  Pc  is  possible,  it  re¬ 
quires  extensive  computation  that  might  be  infeasible 
in  practice  or  not  worthwhile.  Therefore,  finding  easily 
computable  bounds  on  Pc  become  of  practical  impor¬ 
tance. 


3.3  Computing  upper  and  lower  bounds  on  Pc 
*  3,3.1  Upper  bound  using  |detG|:  A  well- 

known  result  of  multi-dimensional  geometry  is  that  the 
volume  of  the  parallelopiped  formed  by  the  basis  vec¬ 
tors  of  a  (nondegenerate)  lattice  is  given  by  |detG|. 
Since  the  Voronoi  cells  of  the  lattice  points  partition 
the  space  with  the  same  density  (cells  per  unit  volume) 
as  the  parallelopiped  cells,  it  follows  that  the  volume  of 
the  Voronoi  cells  are  also  equal  to  |detG|.  Therefore, 
|detG|  which  is  the  volume  of  the  region  for  integrat¬ 
ing  the  noise  probability  density  function,  gives  an  idea 
of  Pc,  and,  roughly  speaking,  the  larger  this  value,  the 
larger  the  probability  of  correct  integer  estimation.  Of 
course,  the  shape  of  the  Voronoi  cell  is  also  a  major  fac¬ 
tor  in  Pc,  however,  when  the  variance  of  noise  in  every 
dimension  is  equal,  an  upper  bound  for  Pc  is  found  if 
we  assume  the  Voronoi  cell  is  a  ^'-dimensional  sphere. 

The  volume  of  a  g-dimensional  sphere  of  radius  p  is 
aqp^  where  Qg  =  7r^/^/r(g/2  +  1)  (r(-)  is  the  Gamma 
function.)  Therefore,  the  radius  of  a  g-dimensional 
sphere  with  volume  |detG|  is  p  =  -^|det  G|/ag  and 
an  upper  bound  on  Pc  becomes 

Pc, up  =  Prob  ^llull  <  ^|det 

with  V  ~  N{0,lqxq)-  The  sum  of  squares  of  q  inde¬ 
pendent  zero-mean,  unit  variance  normally  distributed 
random  variables  is  a  distribution  with  q  degrees 
of  freedom.  If  we  denote  the  cumulative  distribution 
function  (CDF)  of  a  random  variable  with  q  degress 
of  freedom  by  Pj^2  (x^;  q)  we  get 

Fc,up  =  ((|detG|/a,)'/’  ;g)  .  (15) 

This  bound  is  a  very  cheap  one  computationally  since 
it  only  requires  the  determinant  of  G  followed  by  a 
CDF  table  lookup. 

3.3.2  Lower  bound  based  on  dmin:  Given 
the  lattice  L  =  L(G)  with  G  E  the  shortest 

lattice  vector  or  the  minimum  distance  of  the  lattice 
dmin  is  defined  by 

dmin  =  min  ||Gz||.  (16) 

zGZ»,  z^O 

dmin  is  actually  the  distance  between  the  origin  and  its 
closest  neighbor  in  the  lattice.  A  ball  of  radius  dmin/2 
is  the  largest  ball  centered  at  the  origin  that  lies  in  Vo  • 

Clearly,  the  probability  of  noise  falling  in  a  ball  cen¬ 
tered  at  the  origin  and  of  radius  d/2  where  d  <  dmin 
gives  us  a  lower  bound  on  Pc-  This  lower  bound  is  given 
by 

■Pc.iow  =  (^;g^  where  d<  dmin.  (17) 


So  far  we  haven’t  discussed  how  to  compute  dmin-  Un¬ 
fortunately,  computing  dmin  is  conjectured  to  be  NV- 
hard,  and  to  date,  no  polynomial-time  algorithm  has 
been  found  to  compute  dmin-  Therefore,  methods  are 
preferred  that  compute  bounds  on  d^in  that  have  very 
low  computational  complexity.  These  bounds  can  be 
used  in  (17)  to  get  very  fast  bounds  on  Pc-  There  are 
many  ways  to  bound  d^in  as  given  in  [8].  One  method 
is  the  following. 

Gram- Schmidt  method  lower  bound  on  dmin- 
Suppose  G  =  [gi92---9q]  and  is  the 

Gram-Schmidt  orthogonalization  of  (51,52, jff®)  as 
defined  in  (22).  Then 

dmin  >d  =  min(||5l'||,|l5|||,...;||5j||).  (18) 

Re-defining  G  by  finding  a  reduced  basis  for  the  lat¬ 
tice  using  the  LLL  algorithm  would  result  in  a  tighter 
bound  since  the  basis  vectors  (and  hence  the  ^*s)  will 
become  shorter. 


4  Estimation  Problem 

4.1  Nearest  Lattice  Vector  Problem 
The  main  computational  task  of  the  estimation  step  is 
the  solution  to  the  integer  least-squares  problems  (8) 
or  (12).  Problem  (12)  is  known  as  the  nearest  lattice 
vector  problem  and  is  AfV-hard  (cf.  [1].)  However,  in 
practice,  there  are  reasonably  efficient  ways  to  solve 
this  problem  which  is  the  main  topic  of  this  section. 

When  T  >  0  in  (8)  is  diagonal,  can  be  simply  found 
by  rounding  the  components  of  z  to  the  nearest  integer 
since  the  integer  variables  are  seperable  as 

(z  -  z)^T-^(z  -  z)  =  Y^{zi  -  Zif/Tii, 

i=l 

where  Ta  is  the  zth  diagonal  entry  of  T.  A  diagonal 
T  corresponds  to  a  G  with  orthogonal  columns.  In¬ 
tuitively,  this  observation  suggests  that  if  T  is  almost 
diagonal,  or  equivalently,  the  columns  of  G  are  almost 
orthogonal,  rounding  off  z  in  (8)  would  give  a  close  ap¬ 
proximation  to  2;ml  *  However,  G  is  not  always  initially 
given  as  almost  orthogonal  even  if  the  lattice  L  =  L(G) 
is  orthogonal. 

As  noted  before,  the  LLL  algorithm  is  a  useful  tool 
for  finding  an  almost  orthogonal  basis  for  a  given  lat¬ 
tice  L  —  L(G).  Rounding  can  then  be  applied  to  z 
to  get  a  hopefully  good  approximation  for  2:ml-  This 
approximation  can  serve  as  a  good  initial  guess  for  Zml 
in  any  algorithm  for  solving  the  global  optimization 
problems  (8)  or  (12).  A  similar  idea  can  be  found  in  [6] 
and  [7]. 


4.2  Suboptimal  Polynomial-Time  Algorithms 

In  Xhis  subsection,  we  address  suboptimal  polynomial¬ 
time  algorithms  for  calculating  the  minimum  of  \\y  — 
Gz\\  over  the  integer  lattice.  These  suboptimal  algo¬ 
rithms  are  important  for  a  few  reasons.  First,  they  can 
be  performed  efficiently  with  a  guaranteed  low  worst 
case  complexity.  Secondly,  suboptimal  algorithms  pro¬ 
vide  a  relatively  good  initial  guess  or  relatively  tight 
upper  bound  for  any  global  optimization  algorithm. 
Finally,  these  suboptimal  algorithms  might  find  the 
global  optimum  as  they  often  do  in  practice.  If  any 
lower  bound  d  on  dmin  (d  <  dmin)  is  known,  a  sufficient 
condition  for  the  suboptimal  minimizer  Zgub  to  be  the 
global  minimizer  z^h  is  simply  given  by 

\\y  ~  G^^subll  ^  2  - ^  ^sub  =  ^MLj  (19) 

as  there  is  only  one  lattice  point  in  a  ball  centered  at 
G^^sub  and  with  radius  dmin/2.  Note  that  d  is  usu¬ 
ally  a  byproduct  of  the  verification  step,  and  therefore, 
condition  (19)  can  be  checked  without  any  additional 
computation  for  finding  d. 

One  suboptimal  polynomial-time  algorithm  is  the  fol¬ 
lowing  (cf.  [4].) 

Suboptimal  algorithm  for  finding  an  approxi¬ 
mately  nearest  lattice  point  based  on  rounding. 
Suppose  G  6  and  y  G  are  given.  A  subopti¬ 
mal  solution  Zsub  €  in  the  sense  that 

\\y  -  Gzsubll  <  (1  +  25  (4.5)*/^)  min  II5  -  Gz||, 

^  z  e  Z« 

(20) 

exists  and  is  given  by^  Zsub  =  F\G  ^y\  where  G  is 
the  lattice  generator  matrix  after  performing  the  LLL 
algorithm  on  G  and  G  =  GF. 

Note  that  although  the  worst  case  bound  in  (20)  ap¬ 
pears  to  be  loose,  this  suboptimal  algorithm  works 
much  better  in  practice  as  reported  in  the  literature 
(cf.  [5].) 

Another  suboptimal  polynomial-time  algorithm  for 
finding  an  approximately  nearest  lattice  point  is  due 
to  Babai  (1986)  (cf.  [1],  [4]  and  [8].) 

4.3  Searching  for  Integral  Points  Inside  an  El¬ 
lipsoid 

Once  a  candidate  (or  guess)  z  =  z  to  the  minimizer  of 
||y— Gz||  is  found,  we  need  to  check  whether  there  exists 
anyother  z  e  satisfying  ||y  “  Gz||  <  \\y  -  Gz\\.  If  no 
such  z  exists  then  z  is  the  global  minimizer,  otherwise, 
we  can  find  a  better  candidate  for  the  minimum  of 
Gz||  that  we  need  to  check  for  its  global  minimality  as 
well  and  so  on.  Therefore,  an  important  step  in  finding 

^  [-J  is  the  componentwise  rounding  operation  to  the  nearest 
integer. 


the  minimum  of  \\y-Gz\\  is  searching  for  integral  points 
(points  with  integer  coordinates  z)  inside  an  ellipsoid 
||y-G2:||  <  r.  Refer  to  [8]  for  a  method  to  perform  this 
exhaustive  search.  It  is  shown  that  by  putting  G  into  a 
lower  triangular  form  using  a  unitary  transformation, 
this  search  can  be  performed  more  efficiently. 

4.4  Global  Optimization  Algorithm 

Basically,  the  global  optimization  algorithm  for  finding 
the  minimizer  of  \\y  -  Gz\\  consists  of  a  good  initial 
guess  (using  suboptimal  algorithms  of  §4.2)  followed 
by  an  efficient  e^aMStive  search  (§4.3.)  Refer  to  [8]  for 
details. 

4.5  Summary 

In  this  section  we  sketched  a  method  for  solving  the  in¬ 
teger  least-squares  problem  for  resolving  the  integer  pa¬ 
rameters  in  the  linear  model.  In  practice,  this  method 
is  very  efficient  for  problems  with  a  few  ten  integer  vari¬ 
ables.  The  success  of  this  method  mainly  relies  on  a 
guaranteed  reasonably  good  initial  guess  for  the  mini¬ 
mizer  by  using  the  LLL  algorithm. 

After  2  is  estimated  by  the  global  optimization  algo¬ 
rithm  in  §4.4,  the  maximum  likelihood  estimate  of  x, 
the  real  parameter,  is  straightforward.  Xml  can  be 
found  as  in  (9). 


5  Extensions 

In  this  section  we  point  out  extensions  to  model  (1). 
If  there  is  a  priori  knowledge  on  the  distribution  of  x, 
similar  formulas  as  in  §2  can  be  given  for  the  {maximum 
a  posteriori)  estimates  (cf.  [8].)  Now  suppose  that  we 
assume  a  state-space  structure  for  the  real  parameter 
x'^  =  [xj  •  •  -x^]  as 

=  Fk^k  +  GkUk,  C(0)  =  Co, 

Xk  Hk^k  +  Vk  for  k  =  0,1,2,...,  N 

and  the  observation  y^  is  assumed  to  depend  on  the 
unknowns  Xk  €  and  2  G  R^  through  the  linear 
relationship 

yk  =  AkXk  +  JkZ-\-Wk 

in  which  Co,  and  Wk  are  Gaussian  variables.  In 

this  case,  the  matrices  G  and  y  in  (12)  can  be  updated 
recursively  such  that 

=  argmin  -  G^'‘hl  k  =  l,...,N.  (21) 
zeZ‘1 

This  recursive  method  (similar  to  the  Kalman  filter) 
has  many  numerical  advantages  and  reduces  data  stor¬ 
age  for  applications  in  which  the  data  comes  in  se¬ 
quentially  {e.g.,  GPS  surveying.)  A  standard  two-pass 


smoothing  algorithm  (cf.  [3])  can  be  used  to  find  the 
estimates  of  x^  for  i  once  2ml  is  computed 

from  (21).  Refer  to  [8]  for  details. 

6  Simulations 

The  simulations  in  this  section  are  based  on  a  setup 
similar  to  the  global  positioning  system  (GPS)  but  in 
two  dimensions  (Figure  1.)  In  general,  navigation  us¬ 
ing  phase  measurements  from  sinusoidal  signals  trans¬ 
mitted  from  far  away  locations  in  space  with  known 
coordinates  can  always  be  cast  as  an  integer  parameter 
estimation  problem  in  a  linear  model.  The  unknown  in¬ 
tegers  enter  the  equations  as  the  number  of  carrier  sig¬ 
nal  cycles  between  the  receiver  and  the  satellites  when 
the  carrier  signal  is  initially  phase  locked. 


Figure  1:  Simulations  in  this  section  are  based  on  a  simplified 
GPS  setup  in  two  dimensions.  The  radius  of  the  the 
shaded  circle  (earth)  is  6357km.  Three  satellites 
are  assumed  that  orbit  the  earth  at  an  altitude  of 
20200km  and  period  of  12  hours. 

In  our  simulation,  we  have  assumed  a  constant  receiver 
located  at  x^  =  [--50  100]  (of  course  this  is  not  known 
a  priori.)  However,  x  is  known  to  be  Gaussian  with 
mean  zero  and  variance  ax  =  100m  in  both  the  xi  and 
X2  directions.  The  wavelength  of  the  carrier  signal  and 
angular  velocity  for  all  three  satellites  is  A  =  0.19m  and 
Lj  =  l/120®sec~^  respectively.  The  satellites  make  an¬ 
gles  of  90°,  120°  and  45°  with  the  xi  axis  initially  and 
the  direction  of  rotation  for  all  of  them  is  clockwise. 
The  variance  of  phase  measurement  noise  in  units  of 
length  is  taken  as  a  =  0.01m.  The  receiver  measures 
the  (sinusoidal)  carrier  signal  phase  from  each  of  these 
satellites  every  T  =  2  seconds  for  a  period  of  200  sec¬ 
onds.  Refer  to  [2]  and  [8]  for  details. 

Figure  2  gives  the  exact  Pc  (computed  using  Monte 
Carlo  with  1500  random  variates),  and  the  upper  and 
lower  bounds  Pc, up  and  Pc,iow  (lower  bound  using 
Gram-Schmidt  method  lower  bound  on  d^j^)  as  a  func- 


tion  of  time.  Clearly,  the  upper  bound  Pc, up  appears 
to  be  a  very  good  approximation  to  Pc  and  as  Pc  1, 
the  lower  bound  Pc, low  becomes  exact  as  well.  This  is  a 
very  good  feature  as  we  are  not  conservative  in  bound¬ 
ing  Pc  when  Pc  is  high  and  there  is  high  reliability  in 
our  estimates. 


Figure  2:  Fc  {dash-dot  curve),  Pc, up  (upper  bound  using 
|detG|;  upper  solid  curve),  Pc.iow  (lower  bound 
using  lower  bound  on  dmin  using  Gram- Schmidt 
method;  lower  solid  curve)  as  a  function  of  time. 


In  practice  we  might  never  compute  the  exact  Pc  and 
we  have  to  live  with  (the  inexact  but  much  more  easily 
computable)  bounds.  Given  these  bounds  as  in  Fig¬ 
ure  2,  we  can  conclude  that  Pc  is  not  greater  than  0.99 
after  132  seconds,  so  with  a  confidence  level  of  99%,  our 
reliability  on  Zml  is  low.  We  have  to  wait  for  another 
30  seconds  until  the  lower  bound  hits  the  Pc  =  0.99 
line.  Prom  hereon,  it  is  guaranteed  that  Pc  >  0.99  and 
therefore  is  reliable. 

We  have  plot  the  number  of  inner  iterations  in  the 
global  optimization  algorithm  to  compute  Zml  vs.  time 
in  Figure  3.  A  number  of  iterations  equal  to  zero 
means  that  the  initial  value  for  z  found  in  the  algo¬ 
rithm  (by  performing  the  LLL  algorithm  and  rounding 
off)  is  guaranteed  to  be  the  global  minimizer  accord¬ 
ing  to  (19)  where  d  is  the  lower  bound  found  on  dmin 
using  the  Gram-Schmidt  method  in  §3.3.2.  Note  the 
low  number  of  iterations  (cf.  [8]  for  details.)  In  fact, 
the  global  optimization  algorithm  is  very  efficient  in 
practice  and  can  be  solved  instantly  by  a  computer. 

It  is  very  interesting  to  note  that  for  high  Pc  (Pc  >  0.9 
or  after  95  seconds),  the  only  case  for  which  we  are 
really  interested  in  Zml  because  of  high  reliability,  the 
number  of  inner  iterations  is  zero  and  therefore  the 
global  optimization  algorithm  is  extremely  efficient.  In 
general  this  is  true,  and  the  efficiency  of  the  global 
optimization  algorithm  is  greatly  enhanced  as  Pc  1. 


Figure  3:  (a)  Pc  as  a  function  of  time.  (6)  Number  of  in¬ 
ner  iterations  in  the  global  optimization  algorithm 
for  finding  zml  as  a  function  of  time.  For  high  Pc, 
the  initial  guess  for  the  global  minimizer  is  guaran¬ 
teed  to  be  the  global  minimizer  using  a  very  cheaply 
computed  lower  bound  on  dm  in  (using  the  Gram- 
Schmidt  method.)  Therefore,  the  global  optimiza¬ 
tion  algorithm  becomes  extremely  efficient  for  high 
Pc 

Figure  4  shows  Pc  as  a  function  of  time  and  the  times  at 
which  the  integer  parameter  z  is  resolved  correctly  us¬ 
ing  global  optimization  and  the  times  at  which  it  is  re¬ 
solved  correctly  using  rounding  off  (z.e.,  Zest  =  J ; 

without  using  the  LLL  algorithm  of  course.)  The  re¬ 
sults  clearly  show  that  simple  rounding  shouldn’t  be 
used  in  practice  since  even  after  200  seconds  we  are 
still  not  able  to  resolve  z  by  this  method.  On  the  other 
hand,  using  global  optimization,  z  is  resolved  after  90 
seconds  or  when  Pc  >  0.8. 


7  Conclusions 

In  this  paper  we  considered  parameter  estimation  in 
linear  models  when  some  of  the  parameters  are  known 
to  be  integers.  Simulation  results  show  that  if  we  ne¬ 
glect  the  integer  nature  of  the  parameters  (treat  these 
parameters  as  being  real  and  then  round  off),  we  get 
very  inexact  estimates  in  practice.  The  main  compu¬ 
tational  effort  in  the  estimation  step  turns  out  to  be 
the  solution  to  an  integer  least-squares  problem  which 
is  in  fact  A/^P-hard.  Verifying  the  maximum  likelihood 
estimates  seems  to  be  a  problem  as  hard  as  the  esti¬ 
mation  step  (conjectured  to  be  ATP-hard.)  The  prob¬ 
ability  of  correct  integer  parameter  estimation  Pc  was 
chosen  as  a  useful  reliability  measure  for  our  estimates. 
A  polynomial-time  algorithm  due  to  Lenstra,  Lenstra 
and  Lovasz  (LLL  algorithm)  found  to  be  very  useful  in 
the  estimation  and  verification  problems. 


Very  easily  computable  bounds  can  be  found  on  Pc  and 


Figure  4:  Simulation  result  when  =  [0  0].  (a)  Pc  {dash-dot 
curve)  and  the  times  at  which  z  is  resolved  correctly 
using  global  optimization  {solid  curve^  a  high  means 
that  z  is  resolved  correctly  or  zml  =  z,  while  a  low 
means  that  it  is  not)  as  a  function  of  time.  (6)  Pc 
and  the  times  at  which  z  is  resolved  correctly  using 
simple  rounding  as  a  function  of  time. 


these  bounds  become  exact  as  Pc  is  close  to  one;  the 
only  case  that  we  are  really  interested  in  our  param¬ 
eter  estimates  because  of  their  high  reliability.  The 
global  optimization  algorithm  for  solving  the  integer 
least-squares  problem  works  very  well  for  problem  sizes 
of  a  few  tens  of  variables  which  is  typical  in  applica¬ 
tions  such  as  GPS.  Moreover,  the  global  optimization 
algorithm  becomes  even  more  efficient  for  high  Pc- 

These  observations  suggest  the  following  general  out¬ 
line  for  integer  parameter  estimation  in  linear  models. 

General  outline  for  integer  parameter  estima¬ 
tion  in  linear  models.  Under  the  setup  of  §1: 

•  Update  G  after  every  measurement,  say  recur¬ 
sively,  as  noted  in  §5. 

•  Compute  Pc, up  in  (15)  which  is  the  upper  bound 
based  on  the  determinant  of  G.  In  practice,  this 
bound  is  very  close  to  Pc  as  demonstarted  in  the 
simulations. 

•  When  Pc, up  is  large  enough  (say  Pc, up  >  0.99), 
compute  the  lower  bound  Pc, low  in  which  the 
Gram-Schmidt  method  is  used  to  provide  a  lower 
bound  on  dmin  (apply  the  LLL  algorithm  on  G  to 
get  a  tighter  bound.) 

•  When  Pc, low  is  high  (say  Pc,iow  >  0.99),  we  can 
assume  that  z  =  Zml-  Therefore,  only  at  this  step 
eventually  solve  the  global  optimization  prob¬ 
lem  (12).  Since  Pc  is  close  to  unity,  the  algorithm 
for  solving  (12)  would  be  extremely  efficient. 


•  Once  z  =  has  been  resolved,  there  is  no  in- 
^  tegef  parameter  to  worry  about.  Use  standard 
methods  to  estimate  the  real  parameter  x  (e.^., 
two-pass  smoothing  algorithm,) 
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Appendix:  Lenstra,  Lenstra,  Lovasz  (LLL) 
Algorithm 


Suppose  a  lattice  L  =  L(G)  with  G  =  [pi  p2  **  *59]  is  given.  A 
reduced  basis  for  L  can  be  obtained  as  follows, 


•  Perform  the  Gram-Schmidt  orthogonalization  procedure 
on  the  vectors  pi,p2) •  •  •  > **e.,  compute  the  vectors 
9i^92i  •  ‘ '  i9q  through  the  recursion 


ji  —  1  rp 


f*  =  fli  -  5]  1^9.*  fo--  J  = 


i=l 


(22) 


{p*  5  ^2  ’  •  •  ■  ’  ^9  }  orthogonal  basis  for  the  sub¬ 

space  spanned  by  pi  ,p2,  •  •  •  From  the  definition  (22), 
it  is  trivial  that  any  vector  gj  can  be,  expressed  as  a  linear 
combination  of  the  vectors  pj ,  P2  ’  ■  *  *  ’ 


j 

9}  =  ^  .  (23) 

i=l 

where  pji  =  pJp*/||pj|P  for  i  =  1,2, ...  -  1  and  with 

pjj  =  1. 

•  For  j  =  1, 2, . . . ,  p,  and,  given  j,  for  i  =  1, 2, ...  ,j  -  1,  re¬ 
place  Qj  by  gj  -  \pLjii9i,  where  is  the  integer  nearest 
to  pji, 

•  If  there  is  a  subscript  j  violating 

||Pj+i  (^'^) 

then  interchange  gj  and  pj+i  and  return  to  the  first  step, 
otherwise,  stop.  G  =  [pi  p2  ■  •  •  pq]  is  the  reduced  generator 
matrix  of  the  lattice  L. 


The  choice  of  3/4  in  (24)  as  the  allowed  factor  of  decrease  is  arbi¬ 
trary:  any  number  between  1/4  and  1  would  do  just  as  well.  The 
most  natural  choice  (giving  a  best  upper  bound  on  the  product 
Ilpi||llff2||  •  •  •  11^9  II)  would  be  1  instead  of  3/4;  but  the  polynomi¬ 
ality  of  the  resulting  algorithm  cannot  be  guaranteed. 


